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Abstract 

In this work we investigate the quantum dynamics of a model for two single- 
mode Bose-Einstein condensates which are coupled via Josephson tunneling. Using 
direct numerical diagonalisation of the Hamiltonian, we compute the time evolution 
of the expectation value for the relative particle number across a wide range of 
couplings. Our analysis shows that the system exhibits rich and complex behaviours 
varying between harmonic and non-harmonic oscillations, particularly around the 
threshold coupling between the delocalised and self-trapping phases. We show that 
these behaviours are dependent on both the initial state of the system as well as 
regime of the coupling. In addition, a study of the dynamics for the variance of 
the relative particle number expectation and the entanglement for different initial 
states is presented in detail. 
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1 Introduction 



The phenomenon of Bose-Einstein condensation, while predicted long ago [1,2], is nowa- 
days responsible for many of the current perspectives on the potential applications of 
quantum systems. This point of view has arisen with the experimental observation of 
condensation in systems of ultracold dilute alkali gases, realised by several research groups 
using magnetic traps with various kinds of confining geometries. Reviews of the break- 
throughs that have lead to the current state of development can be found in [3,4]. These 
types of experimental apparatus open up the possibility for studying dynamical regimes at 
the frontier between the quantum and classical scenarios, where new macroscopic quantum 
phenomenon can emerge; for example, coherent atomic lasers [5], and the new chemistry 
of atomic- molecular condensates [6]. Subsequently, Bose-Einstein condensates are seen 
as one of the main tools to investigate, verify and improve our understanding of many 
concepts and principles in quantum physics, such as entanglement [7]. It is widely ac- 
cepted that entanglement is the main resource needed for the implementation of quantum 
computation [8]. 

The work on ultracold atomic gases has demonstrated the occurence of Rabi oscilla- 
tions due to quantum tunneling between two internal states of a condensate, similar to 
the Josephson effect in superconductors. Moreover, quite unexpected results were found 
in [9, 10], where it was shown that the temporal evolution of the expectation value of the 
number of particles exhibits a collapse and revival behaviour (see also e.g. [11]). This is 
a novel effect, not observed with spatial and temporal resolution in other systems where 
quantum effects take a macroscopic scale, like the superfiuid and the superconducting 
systems. 

From the theoretical point of view tunneling in Bose-Einstein condensates has been 
widely investigated using a simple two-mode Hamiltonian (see eq. in the next section). 
This model has been studied by many authors using a variety of methods, such as the 
Gross-Pitaevskii approximation [12], mean-field theory [13,14], quantum phase model 
[15], and the Bethe ansatz method [16,17]. Our aim in this work is to expand on the 
theoretic knowledge of tunneling in Bose-Einstein condensates by undertaking a detailed 
and systematic analysis of the quantum dynamics of this two-mode Hamiltonian by means 
of the method of direct numerical diagonalisation. In particular, we will consider the 
temporal evolution of the expectation value for the relative number of particles between 
the two condensates for different choices of the coupling parameters and distinct initial 
states. In this procedure we will keep the total number of particles fixed, while the 
coupling of the Hamiltonian is varied. Using the numerical diagonalisation method, we 
are able to study the model across all couphng regimes. 

We also implement this same approach to investigate the temporal evolution of en- 
tanglement for this model, and contrast this against the variance of the relative number 
of particles. As discussed in [18], there are generally two approaches for creating entan- 
glement in many-body systems. One is to engineer a gapped system whose ground state 
is known to be entangled. By then sufficiently cooling the system to the ground state 
configuration, it will be entangled. In this respect there have been studies of ground 
state entanglement for the XYh model [19,20], the BCS model [21] and two-mode Bose- 
Einstein condensates including the present model [14] . Alternatively, one can manipulate 
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a system so that a given initial state temporally evolves into an entangled state. From this 
point of view entanglement dynamics for coupled Bose-Einstein condensates have been 
investigated in [22], using the same model we study here, but with severe constraints on 
the choice of couplings, and also in [14]. Here we will investigate the behaviour of the 
entanglement evolution across different coupling regimes. 

2 The model 

We will study the dynamics of a model for two single-mode Bose-Einstein condensates 
which are coupled via Josephson tunneling. The model is not only applicable to tunneling 
in atomic Bose-Einstein condensates, but also mesoscopic solid state Josephson junctions 
[12,15,23] and non-linear optics [22]. The Hamiltonian is given by 

H = ^(iVi - N,f - ^{N, - N,) - %{a\a^ + a\a^) (1) 
o Z Z 

where we follow the notational conventions of [12] for the couplings. Above, aj,aj [i = 
1, 2) are the creation and annihilation operators for bosons occupying one of two modes, 
labelled 1 and 2, and A^i and N2 are the respective number operators. The operator N = 
Ni + N2 is the total boson number operator and it is conserved. The coupling k provides 
the strength of the scattering interaction between bosons, A/x is the external potential 
and S J is the coupling for the tunneling. We remark that a change of sign in the tunneling 
coupling £j —£j is equivalent to the unitary transformation CLl — > 0,1, (32 — ^ — CL2- In 
that which follows we will mostly discuss the case A/x = 0. 

Let us now, for future use, distinguish different coupling regimes of the model charac- 
terised by different ratios of k/Sj. We first mention that there exists a threshold coupling 
k/Sj = 4/N which signifies the transition between delocalisation and self-trapping [13] 
(see also [18]): 

• Delocalisation — > k/Sj < A/N 

• Self-trapping — > A/N < k/Sj 

Following [12], it is also useful to consider the following three regimes: 

• Rabi regime — > k/Sj << 1/N 

• Josephson regime — > 1/N « k/£j « N 

• Fock regime — > k/Sj » N. 

For these three regimes it is known that an analogy can be drawn between the dynamics 
of ([Q) and the dynamics of certain pendulums [12]. Comparing the above classifications it 
is seen that the threshold coupling lies in the crossover between the Rabi and Josephson 
regimes, thus offering a candidate to precisely define the boundary between these regimes. 

Below the quantum dynamics of ((T)) will be investigated using the method of direct 
numerical diagonalisation, which allows for a study across all couplings. It is well known 
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that the time evolution of any physical quantity is determined through the temporal 
operator U, given by 



TV 

U = J2e-'^'-'\iJn)M (2) 

n=0 

where {A„} and {|'?/'n)} are the eigenvalues and eigenvectors of the Hamiltonian (Q). There- 
fore, the temporal evolution of any state can be calculated by 

TV 

|V'(t)) = ^a„e-^^"*|V'„), (3) 

n=0 

where a„ = (V'nl^) and |0) represents the initial state. Using eqs. ^ and Q we can 
compute the expectation value and the variance (which is a measure of the quantum 
fluctuations) of any physical quantity represented by a self-adjoint operator, as well as 
the entanglement. In particular, the temporal dependence of the expectation value of any 
operator A can be computed using the expression 

{A) = {m\A\m), 

while the variance is given by 

A{A) = (A^) - {Af. 

Following [24], for any pure state density matrix p of a bipartite system the entangle- 
ment is defined by 

E{p) = -tr(pilog2Pi) = -tr(p2log2P2) 

where pi is the reduced density matrix obtained by taking the partial trace of p over the 
subsystem 2. The definition for p2 is analogous. Here we follow the approach argued 
in [14]. Due to indistinguishability, it is not useful to consider entanglement between 
individual bosons. Rather, we take the two modes of the Hamiltonian (P) to define two 
subsystems, which can be distinguished through the operators Ni and N2 representing 
measurement of each subsystem population. In this case, E{p) is equivalent to [14] 

N 



E{p) = -Y.\^n{t)\'iog,{\c^m 



n=0 



where the Cn{t) are the co-efficients of a general state of the system 



TV 



I^W) = $^c„(t)|iV-n,n). 



n=0 



The states \m,n), which collectively provide a basis for the Fock space, are given by 

\m,n) = ^ L- r- |0)' 
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and |0) is the Fock space vacuum. For a system of total particles the maximal entan- 
glement is Emax = log2(iV +1). 



N = lOO 

N = 400 
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Figure 1: Time evolution of the expectation value for the relative number of particles 
for different ratios of the coupling k/Sj from the top (Rabi regime) to the bottom (Fock 
regime): k/£j = 1/N, 1, iV, iV^ for N = 100, 400 and the initial state is |iV, 0). 



3 Expectation value of the relative particle number 

First we investigate the quantum dynamics of the relative particle number operator Ni — 
N2 (or the imbalance) , in the three different regimes, Rabi, Josephson and Fock, as well 
as in the intermediate ones. For this purpose we use the expressions of the previous 
section together with the eigenvalues and eigenvectors obtained directly from numerical 
diagonalisation of the Hamiltonian (0). 

From fig. ^ it is apparent that the qualitative behaviour in each region, using the 
same initial state, does not depend on the number of particles. (In this figure, as in the 
others following, the chosen time interval is one which most clearly shows the relevant 
dynamical behaviour.) We found that in the interval k/Sj G (close to the 

Rabi regime) the collapse and revival time takes the constant value tcr = 47r^. Also, in the 
region where k/Sj — 1 is small and negative, the expectation value still displays collapse 
and revival behaviour with harmonic oscillations occurring at the value k/Sj = 1. At 

^The ratio k/£j — l/N"^ means that we are using k — 1 and £j = N"^ and similarly for the other 
cases. This convention, which will be used throughout the paper unless noted otherwise, fixes the time 
scale. 
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this coupling the period of oscillations is not independent of A^. Still with reference to fig. 
m we observed in the region where k/£j — 1 is small and positive that the expectation 
value has small amplitude periodic oscillation close to the initial expectation value of 1. 
The amplitude of the oscillations disappear as the ratio k/£j becomes larger (near the 
Fock regime) at which point the bosons remain localised to the subsystem of the initial 
configuration. However this transition from small amplitude harmonic oscillations to 
complete localisation appears to be smooth. Hereafter we will focus most of our attention 
to the Rabi and Josephson regimes, including the crossover, as this is where the most 
complex behaviour is to be found. 




Figure 2: Time evolution for the expectation value of the relative number of particles 
for the initial state |0,A^). The dot-dashed line is for = 100 and the solid line is for 
= 400 with Sj = l{mdk = 8/N. 



In fig. 121 we plot the expectation value for the relative number of particles for two 
different cases: = 100 and A^ = 400, with the initial state |0,A^). Using the relation 
k/Sj = 8/N we reproduce figure 5 of [13], taking into account the different notational 
conventions and using a different time interval to show the collapse and revival for both 
cases. For this case we have followed the conventions used in [13] for the time scale by 
setting £j = 1 and k = 8/N. It is apparent that under this convention the collapse and 
revival time is not independent of A^. However, introducing a rescaled collapse and revival 
time t*^ = 8tcr/N, which corresponds to our adopted time scale convention, shows that 
t*j. is independent of A^. A distinguishing feature of fig. El is that the mean value (i.e. 
time-averaged value) of the imbalance is not zero, in contrast to the collapse and revival 
sequences shown in fig. ^ This is because the respective couplings lie in regimes on 
different sides of the threshold value k/Sj = 4:/N, and these cases respectively illustrate 
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the behaviour associated with self-trapping and delocahsation. 

In fig. El we can see in detail for the case = 100 the evolution of the dynamics 
from a collapse and revival sequence with tcr = 47r for k/Sj < 4:/N, through the self- 
trapping transition at k/Sj = A/N, and toward small amplitude harmonic oscillations in 
the imbalance of the localised state when k/Sj = 1. Where there is clear collapse and 
revival in the self-trapping phase we find that tcr increases with increasing k/Sj. Further 
increases in k/£j lead to a decaying of the collapse and revival sequence toward harmonic 
oscillations which occur at k/Sj = 1. 




t t 

Figure 3: Time evolution of the expectation value between k/Sj = 1/N and k/Sj = 1. 
On the left, from top to bottom k/Sj = 1/N,2/N,3/N,4:/N and on the right, from top 
to bottom k/£j = 5/N, 10/N, 50/N, 1, where N = 100 and the initial state is |A^, 0). 



Next we turn our attention to study the evolution of the expectation value {N1—N2) /N 
using a range of different initial conditions. In the first of these studies fig. |3J we illustrate 
the types of different behaviours that occur by using different initial states. Here, it is 
clear that the collapse and revival time depends on the choice of initial state. As expected 
the mean amplitude of the oscillations is related to the imbalance of the initial state: the 
mean amplitude is maximal when the initial state has all particles in the same subsystem 
(|A^, 0) and |0,A^)), and it is minimal when the initial state has the same number of 
particles in each subsystem {\N/2, N/2)). In both regimes shown in fig. |3]the expectation 
value of the relative number of particles is symmetric about zero, with oscillations in the 
relative number disappearing when the imbalance of the initial state is zero. 

As we vary the coupling ratio past the self-trapping threshold value we find an entirely 
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(a) (b) 




Figure 4: Relative number expectation values for (a) the Rabi regime {k/£j = l/N'^); 
(b) between the Rabi and Josephson regimes {k/Sj = 1/N). The initial states, from top 
to bottom are |100, 0), |90, 10), |74, 26), |60, 40) and |50,50). 



different scenario. In fig. El the expectation value for the relative number of particles 
shows a variety of different behaviours. The dynamics is harmonic when the initial state 
corresponds to all particles in the same subsystem, and the oscillation amplitude is small. 
As the imbalance in the initial state decreases, we first see harmonic modulation of the 
amplitude of the oscillations, in stark contrast to the cases shown in fig. |3 Note here 
that the period of the modulation does not depend on the initial state. It is important to 
mention that in the interval k/Sj G [1/A^, 1], which is not shown, we find that the crossover 
from collapse and revival sequences of oscillations with zero mean to oscillations with mean 
approximately equal to the initial imbalance occurs at k/£j = 4:/N, as expected. 

When using other initial states, we can classify them into two categories: symmetric, 
such as a cat-state 

|0CAT) = ^(|iV,O) + |O,iV)) (4) 
and a maximally entangled state 

1 ^ 

|0me) = r^-— J] |iV - n, n), (5) 

or non-symmetric, such as |A^, 0), |A^ — 10, 10), etc. The expectation value for the number 
of particles of any symmetric initial state is zero, due to symmetry of the Hamiltonian 
under interchange of the two subsystems. In such a case some insight into the behaviour 
of the system can be gleaned by investigating the variance of the imbalance, which will 
be explored in the next section. 



8 



(a) 



(b) 



0.999 
0.81 



0.99999995 







^ lU 13 




5 lO 15 20 

1 ' 1 ' 1 ' 

U/vVWVaaaa/WWxaatvWVv-^ 
1.1,1. 


1 


5 lO 15 20 

1 ' 1 ' 1 ' 


1.1,1. 




0.2 
0. 19992 



0. 1 
t 



Figure 5: Relative number expectation values for (a) the Josephson regime [k/Sj = 1); 
(b) between the Josephson and Fock regimes {k/Sj = N). The initial states are, from 
top to bottom: |100,0), |90,10), |74,26), |60,40) and |50,50). 



When a non-zero external potential A/i is added to the system, some new features 
appear in the various regimes. In the extreme Rabi regime, the imbalance population 
oscillates around zero, similar to the A/x = case, but here the dependence on the 
initial state is different and in particular the oscillation amplitude is not zero using the 
initial state \N/2, N/2). Another important point is that the initial states \ni,n2) and 
\n2,ni) produce the same expectation value behaviour. In the crossover between Rabi 
and Josephson regimes, a different behaviour is found: the oscillation amplitude of the 
number of particles also depends on the initial state, but here the initial states \ni,n2) and 
\n2,ni) do not produce the same expectation value behaviour. The symmetry breaking 
of the Hamiltonian due to the term A/i is prominent in this case. 

In the Josephson and Fock regimes the introduction of the external potential only 
produces a significant difference in comparison to the Afi = case when the initial 
imbalance is small. If we consider an initial state where the imbalance population is 
large, then the two initial states \m,n) and \n,m) produce essentially the same results, 
while when \m — n\ is close to zero the same initial states produce different results. This 
aspect can be understood when we write down the first two terms of the Hamiltonian (^J) 
in the following way: 




(iVi - iVj) 




When the fluctuations in the imbalance are small, we see from the terms in the brackets 
that the external potential will be dominant when 
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A^> -(iVi-iVa). 

This means that the largest effect of the external potential will generally occur when the 
initial state is symmetric. 

4 Imbalance fluctuation and entanglement 

Here we will study the extent to which the evolution of the imbalance variance and the 
entanglement are correlated, and how the choice of initial state, as well as the coupling 
regime, determine the main features of the evolution. On one hand, any state of the 
system (with fixed N) which is unentangled must have zero fluctuation in the imbalance. 
The converse however is not true, as the cat-state (jU has maximal imbalance variance 
but relatively small entanglement, and the maximally entangled state (0) has smaller 
imbalance variance than @. 
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Figure 6: Time evolution of the relative particle number fluctuations (on the left) and 
entanglement (on the right) for the initial state |A^, 0) where N = 100, 400. From top to 
bottom the ratio of couplings used are k/8j = l/N^A/N, W/N, 1. 



For the results shown in fig. ^ we have chosen |A^, 0) as the initial state, which is 
unentangled. In the coupling regime k/Sj G [1/iV^, '^/N] where the tunneling is strong, 
we find that the behaviour of the variance of the relative particle number is similar to 
the collapse and revival character of the expectation value (cf. fig. Q). We can observe 
that the entanglement quickly approaches the maximum value and for all later times the 
system is mostly entangled. At the threshold coupling k/Sj = 4:/N, there is no collapse 
and revival behaviour in the imbalance variance and we see that the mean value starts to 
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decrease. Beyond the threshold couphng, the mean value of both the imbalance variance 
and the entanglement decrease. When the ratio k/£j = 1 is reached, we find small 
oscillations close to zero for both the imbalance variance and the entanglement. From 
here to the Fock regime (not shown), the imbalance variance practically goes to zero and 
the entanglement essentially disappears. 
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Figure 7: Time evolution of the relative particle number fluctuations (on the left) and 
entanglement (on the right) for the initial state |0cat) where = 100, 400. From top to 
bottom the ratio of couplings used are k/Sj = 1/N'^,1/N,A/N,10/N,1. 



For flg. [7| the initial state is the cat-state (@)) which has entanglement E{p) = 1. In 
contrast to flg. El the mean value of the imbalance variance increases beyond the threshold 
coupling. As the ratio k/Sj approaches unity, small amplitude periodic oscillations are 
observed very near to the maximal value of A^^. For the entanglement the picture is 
similar to flg. IHl For k/Sj = 1 the entanglement is found to display periodic oscillations 
about a mean value close to the value for the entanglement of the initial state, which is 
non-zero in this instance. 

Finally, flg. |H1 shows the results obtained when the initial state is the maximally 
entangled state (0). Here, a much different behaviour is found in comparison to the 
previous examples. In the interval k/Sj E [1/A^^, 1/A^] the mean values for the variance 
of the imbalance are smaller than in the previous two examples. At the threshold coupling 
the mean is substantially larger than the mean value below threshold coupling, and the 
mean value again decreases as k/Sj approaches unity. The most striking feature at 
k/Sj = 1 is that there is no discernable periodicity in the imbalance variance, whereas in 
the analogous cases shown in flgs. IHl El there is clear periodicity. For the entanglement, 
it appears that the evolution is not substantially dependent on the coupling regime when 
the initial state is maximally entangled. All cases shown have much the same mean 
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Figure 8: Time evolution of the relative particle number fluctuations (on the left) and 
entanglement (on the right) for the initial state I^me) where N = 100, 400. From top to 
bottom the ratio of couphngs used are k/Sj = l/N'^, l/N^A/N, W/N, 1. 



entanglement over time. Like the imbalance variance, there is also no periodicity in the 
entanglement evolution for the coupling k/Sj = 1. 



5 Conclusion 

To summarise, we have investigated the quantum dynamics of a model for two Josephson- 
coupled Bose-Einstein condensates across a wide range of coupling regimes and various 
initial states, using the method of direct numerical diagonalisation. Our analysis shows 
that despite the apparent simplicity of the Hamiltonian, diverse features including collapse 
and revival of oscillations, amplitude modulated oscillations and non-harmonic behaviour 
are displayed. The most striking changes in the dynamics occur as the system crosses the 
threshold coupling k/£j = A/N. 
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